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ABSTRACT 


The  problem  of  predicting,  on  the  basis  of  an  observed  sample  , 
from  an  inverse  Gaussian  distribution,  the  mean  of  a  future  random 
sample  (or  a  single  future  observation)  from  the  same  distribution  is 
considered.  Approximate  prediction  intervals  are  proposed,  and  their 
accuracy  is  investigated  via  extensive  Monte  Carlo  simulations.  The 
results  are  useful  for  predicting  the  next  first  passage  time  for  a 
Brownian  motion  with  positive  drift  or  the  failure  time  of  an  item 
having  inverse  Gaussian  life  distribution. 
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Key  Words:  Prediction;  First  passage  time  distribution;  Life  testing; 
~  Reliability;  Monte  Carlo  simulations. 


1.  INTRODUCTION 


The  inverse  Gaussian  distribution  arises  as  the  first  passage  tine 
distribution  of  a  Brownian  notion  process  with  drift  (Cox  and  Miller,  1965). 
Tweedie  (1957a,  1957b)  discussed  the  statistical  properties  of  this  distri¬ 
bution  and  noted  the  similarities  between  the  sampling  distribution  theory 
of  the  maximum  likelihood  estimators  of  the  inverse  Gaussian  parameters  and 
those  of  the  mean  and  variance  of  the  normal  distribution.  Recently,  Chhikara 
and  Folks  (1977)  proposed  the  inverse  Gaussian  distribution  as  a  lifetime 
model  in  reliability  studies  which  has  several  advantages  over  the  lognormal 
model  when  there  is  a  high  occurrence  of  early  failures.  Further  distribution 
theory  for  the  inverse  Gaussian  was  given  by  Chhikara  and  Folks  (1974,  1975), 
and  its  use  in  tests  for  positive  drift  in  Brownian  motion  was  discussed  by 
Seshadri  and  Shuster  (1974)  and  Chhikara  and  Folks  (1976),  among  others. 

The  form  of  the  probability  density  function  of  the  two-parameter  inverse 
Gaussian  distribution  generally  used  is 

f (x;y,X)  -  U/(2ttx3)3Js  exp[-X(x-y)2/(2xy2) 3,  x>0,(X,y>0)  .  (1.1) 

Other  parametric  forms  are  summarized  by  Johnson  and  Kotz(1970).  In  (1.1), 

3 

y  is  the  mean  and  X  is  a  shape  parameter.  The  variance  is  y  /X,  so  that  y 
is  not  a  location  parameter  in  the  usual  sense. 

For  the  inverse  Gaussian  lifetime  model,  Fadgett  (1979)  obtained 
approximate  confidence  bounds  for  the  reliability  function,  and  Padgett  and 
Wei  (1979)  added  a  threshold  parameter  or  "guarantee  time"  and  studied  the 
estimation  problem  for  the  resulting  three-parameter  model. 

In  this  paper  a  prediction  Interval  is  proposed  for  the  mean  of  m 
Independent  future  observations  from  (1.1),  Y-,...,Y  ,  on  the  basis  of  a  past 


iniiinfiifriiill 


random  sample  of  size  n,X, . X  ,  from  the  same  distribution.  That  is,  it 

i  n 

is  desired  to  obtain  functions  of  X, »...,X  ,  L(X, ,...,X  )  and  U(X, ,...,X  ), 

l  n  l  n  in 

so  that  for  a  specified  value  Y, 


P[L(X.,...,X)  S  Y  <  U(X. . X  -  1-Y, 

i  n  m  i  n 


(1.2) 


where  Y  ■  J  Y  /m.  If  m*l,  (1.2)  provides  a  prediction  interval  for  a 
®  i-1  1 

single  future  observation  Y.  Note  that  the  correct  interpretation  of  the 

probability  statement  (1.2)  is  that  if  repeated  past  samples  are  taken  and 

used  to  predict  the  mean  of  the  future  samples  by  the  intervals  from 

L(X, ,...,X  )  to  U(X, ,...,X  ),  then  Y  will  be  contained  in  the  intervals  for 
in  in  m 

a  proportion  1-Y  of  such  pairs  of  samples.  The  intervals  defined  by  (1.2) 
will  be  useful  in  predicting  the  next  first  passage  time  of  a  Brownian  motion 
with  positive  drift  or  the  future  failure  time  of  an  item  having  inverse 
Gaussian  life  distribution  based  on  a  past  sample,  for  example. 

Prediction  intervals  for  the  normal  and  exponential  distributions  have 
been  studied  extensively.  For  example,  among  others,  Hewett  (1968)  and 
Lawless  (1970,  1972)  considered  prediction  intervals  for  the  exponential  dis¬ 
tribution,  and  Chew  (1968),  Hahn  (1969,  1970),  Hall  and  Prairie  (1973),  and 
Mann  and  Fertig  (1977)  studied  prediction  intervals  for  normal  distributions. 

It  is  very  difficult,  if  not  impossible,  to  obtain  exact  prediction 
intervals  of  the  form  (1.2)  for  the  inverse  Gaussian  distribution  (1.1)  when 
both  parameters  are  unknown.  In  Section  2,  approximate  prediction  Intervals 

for  Y  are  proposed  based  on  estimates  of  y.  The  performance  of  the  Intervals 

n 

in  terms  of  coverage  probabilities  is  investigated  via  extensive  Monte  Carlo 
simulations,  and  the  results  are  reported  in  Section  3.  The  proposed  approxi¬ 
mate  prediction  intervals  seem  to  perform  very  veil.  An  example  is  given 
in  Section  4. 


*  *  v*-*  i  j *  j  ,j 
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2.  THE  PROPOSED  PREDICTION  INTERVALS 


Let  X, be  a  random  sample  of  size  n  from  the  inverse  Gaussian 
x  n 

distribution  (1.1)  and  let  Y,  ,...,Y  be  a  future  random  sample  of  size  m 

x  m 

from  the  same  distribution,  independent  of  the  first  sample.  The  maximum 

~  n 

likelihood  estimators  of  p  and  X  are  given  by  p  *  X  ■  I  X./n  and 
A_.  n  n  i-l  1 

X  ■  £  (1/X  -1/X  )/n,  respectively  (Tweedie,  1957a,  1957b).  Also,  Tweedle 

i-1  n  _ 

(1957a)  showed  that  X  and  X  are  stochastically  independent,  X  has  inverse 

n  n 

a 

Gaussian  distribution  with  parameters  p  and  nX,  and  nX/X  has  chi-square 
distribution  with  n-1  degrees  of  freedom.  Also,  Y  has  inverse  Gaussian 

B 

distribution  with  parameters  p  and  mX,  and  by  a  result  of  Shuster  (1968), 

mX(Y  -p)2/p2Y_  has  chi-square  distribution  with  one  degree  of  freedom  and 
n  n 

A  ^  _  2  2- 

is  independent  of  nX/X.  Thus,  the  random  variable  (n-l)mX(Y  -p)  / (np  Y  ) 

B  B 

has  F  distribution  with  (l,n-l)  degrees  of  freedom. 

Let  F^(l,n-1)  denote  the  value  such  that  an  F  random  variable  with 
(l,n-l)  degrees  of  freedom  satisfies  P[FSFy(l,n-l)3  ■  1-Y.  Then  from  above 


(Y  -p)  AF  (l,n-l) 

PC— f - *  — * - - - 3  -  1-Y. 

pi  m(n-l)X 

ID 


(2.1) 


If  the  parameter  P  is  known,  then  (2.1)  can  be  solved  for  an  exact  (1-Y) 

prediction  interval  for  Y  .  However,  if  p  is  unknown,  then  an  exact  prediction 

B 

interval  cannot  be  obtained  from  (2.1).  An  approximate  1-Y  prediction 

interval  may  be  obtained  as  follows.  Since  we  have  the  random  sample 

2 

Xj,...,Xb  at  hand,  suppose  that  p  in  the  denominator  in  (2.1)  is  approxi- 
•2 

mated  by  X.,  so  that 
n 


(Y  -p)2  nX2  F  (l,n-l) 

P[  J5 - s  — p  T  - - 

Y_  m(n-l)X 


3  *  1-Y. 


(2.2) 
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Also,  since  y  esn  be  estimated  by  the  future  sample  mean, 

nX  -HnY 


of  vi  in  (2.2)  can  be  approximated  by  — — — 

n+» 


Then 


(Y  -vj)4 

ZD 


x2 


~  (~5r)2  (Y  4  -2-  -  2X). 
n*m  m  ^ 


the  unknown  value 


Hence,  (2.2)  becomes 


X2  (n+m)2X2  F  (l,n-l) 
P[Y  4  —  s  - 2— X - 

"  V  na(n-l)\ 


4  2X  )  m  1-y. 

U 


(2.3) 


The  inequality  in  (2.3)  can  be  solved  for  Y  by  finding  the  roots  of  the 

D 

quadratic  equation 


Y2  -  c(X)Y  +  X2  *  0. 

VI  VI 


(2.4) 


where  c(X)  • 


(n-hn)2x2  F  (ltn-l) 
- 5 _ 1 - - - 


nm 


(n-l)X 


4  2X  .  It  is  easy  to  show  that  both  roots 
n 


of  (2.4)  are  always  real  and  positive  and  that  the  quantity  Y  4  - — 

■  Y 


satisfies  the  inequality  in  (2.3)  when  f  is  between  these  two  roots.  Denote 

0 

the  smaller  root  of  (2.4)  by  «nd  the  larger  root  by  U^QC)  where 

X  -  (X. ,...,X  ).  Then  the  Interval  (L. (X),  0, (X))  provides  an  approximate 
“in  l  —  l  — 

(1-Y)  prediction  interval  for  5  .  The  closeness  of  the  approximation  will 

0 

be  Investigated  in  Section  3  by  Monte  Carlo  simulations.  For  large  n,  since 

X_  is  very  near  y  with  probability  one  by  the  strong  law  of  large  numbers, 
n 

the  prediction  probability  will  be  close  to  the  nominal  value  (1-Y). 

A  second  approximate  prediction  Interval  for  Y_  may  be  obtained  by  using 

0 

•  2  2- 

the  first  few  terms  of  the  Taylor  series  expansion  of  g(y)  •  (Y  -y)  /(y  Y  ) 

0  ffl 

about  X_.  Zn  the  Monte  Carlo  study  reported  in  Section  3,  the  approximation 

U 


^•yr«fv 
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using  only  one  or  two  terms  in  the  series  was  not  very  good.  However,  the 
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approximation  using  three  terms  performs  fairly  well  when  the  variance  y  /X 
of  the  inverse  Gaussian  distribution  is  small  and  n  is  relatively  large. 

The  three  terms  are 


-  X_2[(6-8y/X  +  3y2/X2)  Y  +  X2  /  Y 
n  n  n  m  n  m 

-  6X  +  6y  -  2y2/X  3.  (2.5) 

n  n 

From  (2.1)  and  (2.5), 

P[ (6-8y/X  +  3y2/X2)Y  +  X2/? 

n  n  m  n  m 

nF(l,n-l)X2  ,  _ 

s  — * - - S-  +  6X  -6y  +  2 y4/X  ] 

m(n-l)X  n  n 

~  1-Y.  (2.6) 

The  probability  (2.6)  still  depends  upon  y,  and  as  before,  an  approximation 

must  be  used.  A  lower  confidence  bound  on  y  was  stated  In  Padgett  (1979)*  Using 

Xn  for  y  on  the  left  side  of  the  Inequality  in  (2.6)  and  the  lower  confidence 

bound  [1/X  +  (F  (l,n-l)/((n-l)ji  X))]”2  on  the  right-hand  side  gives  a  rough 
n  y  n 

approximation 

X* 

PC?  +  —  s  c,(X)3  SB  1-Y, 
ms  *  " 


(2.7) 


where 


with 


nFv(l,n-l)X* 

c,(X)  -  —5 l - — -S-  +  2XQ(X), 

1  m(n-l)X 

Q(X)  -  3tX  Fv(l1n-l)]‘5/[((n-l))i)!j  +  (X  Fjl.n-l))1*] 
—  n  y  n  Y 

+  (n-l)X/[  ((n-DX)*5  +  (XnF  (1,b-1))V. 


As  before,  (2.7)  has  a  unique  solution  of  the  form 
P[L_(X)  S  Y  &  U,(X)3  s  1-Y,  where  L-(X)and  U,(X)are  the  real  positive  roots 
of  the  equation  Y^  -  CjMY  +  jP  *  0. 

Other  approximations  for  prediction  intervals  for  Y  are  possible,  but 

n 

the  two  given  here  seem  to  behave  better  than  any  of  the  many  others  tried 
in  the  Monte  Carlo  study.  The  Monte  Carlo  simulation  study  is  described  and 
some  of  the  results  are  reported  in  the  next  section. 


3.  MONTE  CARLO  SIMULATIONS 


Before  describing  the  Monte  Carlo  simulation  results,  a  brief  description 
of  the  procedure  for  generating  a  random  number  from  the  inverse  Gaussian 
distribution  is  given.  Let  X  denote  a  random  variable  with  inverse  Gaussian 
distribution  (1.1).  Then  as  before  Y  "  X(X-y)^/(Xy^)  has  chi-square  dis¬ 
tribution  with  one  degree  of  freedom.  Hence,  a  value  y  of  Y  is  generated, 
and  the  solution  for  x  in  terms  of  y  from  above  is 

x  ■  [y(2X  +  yy)  +  y(yy(yw  +  ^X))**)^.  (3.1) 

It  is  obvious  that  both  solutions  given  by  (3.1)  are  positive  and  that  the 
plus  sign  gives  x>y  and  the  minus  sign  gives  x<y,  with  x»y  only  if  y-0. 


Therefore,  a  value  u  of  a  uniform  random  variable  U  on  (0,1)  is  generated. 
Let  denote  the  solution  from  (3.1)  given  by 


[y(2X  ♦  yy)  +  y(yy(yy  +  4X))*S]/2X 


and  let 


x2  -  Cy(2X  +  yy)  -  y(yy(yy  +  4X))Js]/2X, 


Let  p(y)  *  —  .  Then  choose  the  value  of  the  inverse  Gaussian  random 

variable  X  to  be  x^  if  u  s  p(y)  and  x2  otherwise,  since  the  pdf  of  X  is  given 
by  the  mixture  f^(x)p(y)  +  f2(x)[l-p(y) )  for  a  given  y,  where  f.(x.)  denotes 
the  pdf  of  X1  and  f2(x2)is  the  pdf  of  X2»  Repeating  this  procedure  for  n 
independent  values  of  Y  and  U  yields  a  (pseudo)  random  sample  of  size  n 
from  (1.1). 


In  order  to  investigate  the  behavior  of  the  approximate  prediction 
intervals  proposed  in  Section  2,  extensive  Monte  Carlo  simulations  were 
performed.  For  various  fixed  values  of  Y,  n,  m,  y,  and  X,  estimates  of  the 
coverage  probabilities  and  mean  widths  of  the  prediction  intervals  for  Y 

n 

discussed  in  Section  2  were  computed  based  on  1000  pairs  of  samples  X^,...,XQ 

and  Y, , . . . ,Y  .  The  results  were  essentially  the  same  when  2000  pairs  of 
x  m 

samples  were  used,  so  most  of  the  simulations  were  performed  with  1000  samples 

in  order  to  reduce  computing  time.  For  each  fixed  set  of  values  Y»  n,  m,  X, 

and  y,  1000  pairs  of  samples  x. ,...,x  and  y1t...,y  were  generated  from  the 

Inverse  Gaussian  distribution  (1.1).  For  each  such  pair,  the  approximate 

(1-Y)  prediction  interval,  say  L^(x)  to  U^(x),  was  computed  and  checked  to 

see  whether  or  not  y  was  contained  in  the  computed  Interval.  The  width  of 

n 

the  interval  was  also  computed  and  stored.  Then  the  proportion  of  the  1000 
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such  pairs  of  samples  which  gave  prediction  intervals  containing  the 
corresponding  y^  was  used  as  the  estimate  of  the  actual  coverage  probability. 
Also  the  average  interval  width  was  computed  for  the  1000  pairs  of  samples. 

Tables  1-3  give  some  of  the  Monte  Carlo  results  for  1-Y  *  0.99,  0.95, 
and  0.90  and  various  sample  sizes  n  and  m  and  values  of  y  and  X.  The  symbol 
^  refers  to  the  prediction  interval  (L^X),  1^00)  in  Section  2  and  I2 
refers  to  (l^OO,  U2(X)).  Based  on  the  simulation  results,  1^  has  a  better 
overall  performance  as  an  approximate  (1-Y)  prediction  interval  for  Y  than 

ID 

any  other  approximation  which  was  tried.  The  prediction  interval  12  is 

somewhat  conservative  when  the  variance  of  the  underlying  inverse  Gaussian 
2  3 

distribution,  o  *  y  /X,  is  small  and  does  not  give  a  good  approximation  when 
the  variance  is  large  and  n  is  small,  for  example  y*3,  X«5t,  and  n*5.  As 
expected,  the  approximation  becomes  closer  and  the  average  interval  width 
decreases  as  n  increases.  Also,  as  the  variance  a  *  y  /X  decreases,  the 
prediction  improves  and  the  prediction  intervals  tend  to  be  conservative. 

Both  1^  and  I2  perform  better  for  m»l  than  for  other  values  of  m. 

4.  EXAMPLE 


Chhikara  and  Folks  (1977)  showed  that  the  maintenance  data  reported 
by  Von  Alven  (1964)  on  active  repair  times  (hours)  for  an  airborne  commu¬ 
nications  transceiver  fit  an  Inverse  Gaussian  distribution.  The  n*46  observed 
repair  times  were: 


.2, .3, .5, .5, .5, .5,. 6, .6, .7, .7, .7, .8, .8, 1.0, 1.0, 1.0, 1.0, 1.1, 
1.3, 1.5, 1.5, 1.5, 1.5, 2. 0,2. 0,2. 2, 2. 5, 2. 7, 3. 0,3. 0,3. 3, 3. 3, 


4.0,4.0,4.5,4.7,5.0,5.4,5.4,7.0,7.5,8.8,9.0,10.3,22.0,24.5. 
The  respective  maximum  likelihood  estimates  of  y  and  X”*  were  X^  - 


3.61  and 
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Table  1.  Simulation  Results  for  l-Y  ■  0.99 


V  A 


1  0.25 


1  0.25 


3  0.25 


3  0.25 


1  1 


1  4 


1  4 


1  4 


5  1 


5  4 


5  4 


n  m 


5  1 

5  5 

5  10 


30  1 

30  15 

5  1 

5  5 

5  10 

30  1 

30  15 

15  1 

15  10 

15  20 

5  1 

5  5 

15  1 

15  10 

30  1 

30  15 

15  1 

15  20 

5  1 

5  5 

30  1 

30  15 


Ave .Width 


208.308 

118.512 

135.319 

38.365 

6.640 

3430.014 

1629.692 

2285.851 

368.060 

61.311 

12.354 

4.114 

3.976 

9.157 

5.886 

4.084 

1.671 

3.531 

1.120 

332.821 

94.316 

220.023 

138.620 

62.778 

14.571 


1 

Co  v.Prob. 


0.983 

0.995 

0.990 

0.989 

0.996 

0.977 

0.967 

0.962 

0.990 

0.981 

0.993 

0.994 

0.999 

0.988 

0.997 

0.994 

0.998 

0.991 

0.998 

0.990 

0.992 

0.989 

0.996 

0.992 

0.994 


Ave. Width 


148.0918 

33.834 

19.600 

37.566 

5.313 

2393.380 

420.363 

268.782 

351.081 

36.551 

12.304 

3.641 

3.063 

8.266 

3.993 

4.520 

2.114 

3.933 

1.690 

303.860 

34.503 

166.451 

52.584 

64.342 

15.493 


Co v.Prob. 


0.980 

0.967 

0.940 

0.989 

0.995 

0.967 

0.920 

0.891 

0.990 

0.968 

0.994 

0.994 

1.000 

0.995 

1.000 

0.998 

1.000 

0.996 

1.000 

0.988 

0.970 

0.987 

0.990 

0.994 

1.000 


w 
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Table 

2. 

Simulation  Results 

for  1-T  - 

0.95 

X1 

X2 

p 

\ 

n 

m 

Ave. Width 

Cov.Prob. 

Ave. Width 

Cov . Prob . 

1C 

5 

1 

73.393 

0.942 

53.844 

0.933 

1 

U.  2d 

5 

5 

39.481 

0.952 

13.439 

0.921 

A  A  f* 

30 

1 

21.850 

0.962 

21.815 

0.963 

1 

0.23 

30 

15 

4.421 

0.974 

4.053 

0.976 

a  a  r 

5 

1 

1441.764 

0.948 

1012.451 

0.935 

3 

0.25 

5 

5 

503.665 

0.925 

137.920 

0.858 

3 

0.25 

30 

1 

226.021 

0.957 

217.445 

0.956 

15 

1 

7.161 

0.956 

7.466 

0.966 

1 

1 

15 

10 

2.653 

0.983 

2.787 

0.989 

*1 

50 

1 

5.964 

0.960 

6.359 

0.966 

1 

1 

50 

40 

1.194 

0.990 

1.724 

1.000 

5 

1 

4.420 

0.959 

4.564 

0.981 

1 

4 

5 

5 

2.888 

0.974 

2.704 

0.997 

30 

1 

2.387 

0.944 

2.757 

0.964 

1 

4 

30 

15 

0.810 

0.983 

1.388 

1.000 

15 

1 

198.000 

0.973 

184.395 

0.966 

5 

1 

15 

10 

48.958 

0.965 

31.178 

0.953 

30 

1 

138.490 

0.948 

136.992 

0.950 

5 

1 

30 

15 

24.509 

0.969 

21.461 

0.967 

5 

1 

100.599 

0.956 

81.492 

0.959 

5 

4 

5 

5 

56.788 

0.973 

28.971 

0.974 

30 

1 

37.529 

0.942 

39.581 

0.955 

5 

4 

30 

15 

9.941 

0.980 

12.077 

0.997 

i mutai, 
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Table  3.  Simulation  Results  for  1-Y  *  0.90 


u 

X 

n 

s> 

Ave. Width  x 

Cov.Prob. 

Ave. Width 

2 

Cov.Prob. 

5 

1 

46.122 

0.897 

34.753 

0.888 

1 

0.25 

5 

5 

26.807 

0.931 

10.095 

0.894 

5 

10 

30.468 

0.921 

7.079 

0.868 

5 

1 

739.880 

0.895 

523.699 

0.875 

3 

0.25 

5 

5 

353.523 

0.880 

99.665 

0.815 

5 

10 

494.554 

0.885 

67.965 

0.789 

1 

4 

15 

1 

2.027 

0.893 

2.461 

0.949 

15 

10 

0.948 

0.966 

1.518 

0.999 

5 

4 

15 

1 

31.853 

0.887 

33.583 

0.905 

15 

10 

12.130 

0.950 

13.221 

0.975 

12 


A  .  A 

X  »  0.587  (or  X  ■  1.70).  For  this  data  and  m*l,  a  95%  prediction  interval 
for  the  next  repair  time  is  obtained  from  (1^00,  U^Op)  in  Section  2  as 
(0.3185,  AO. 8393).  Similarly,  (L^QO,  UjOO)  yields  a  95%  prediction  interval 
of  (0.3074,  42.3139),  slightly  wider  than  (1^00,  U^X)).  For  m-10,  a  95% 
prediction  interval  for  the  mean  of  the  next  ten  repair  times  is  found  from 
(1^00,  U2(X))  to  be  (1.2009,  10.8311)  and  (L2(X),  U2(X))  yields  (1.0575,12.2997). 

5.  CONCLUSION 

Prediction  Intervals  for  the  mean  Y  of  m  future  observations  (or  a 

m 

single  future  observation)  based  on  a  current  independent  random  sample  of 

size  n  from  the  same  inverse  Gaussian  distribution  (1.1)  have  been  investigated. 

When  the  mean  y  of  the  inverse  Gaussian  distribution  is  known,  exact  (1-Y) 

prediction  intervals  are  easily  obtained.  If  y  is  unknown,  some  approximate 

(1-Y)  prediction  intervals  for  Y  have  been  proposed.  Based  on  the  results 

m 

of  Monte  Carlo  simulations,  the  approximate  prediction  interval  12  given  by 
(Lj(X),  U^OO )  in  Section  2  is  relatively  simple  to  compute  and  performs 
best  overall  of  those  considered  in  the  study. 
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